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ABSTRACT 

Context. Cosmic rays (CRs) are transported out of the galaxy by diffusion and advection due to streaming along magnetic field 

lines and resonant scattering off self-excited MHD waves. Thus momentum is transferred to the plasma via the frozen-in waves as a 

mediator assisting the thermal pressure in driving a galactic wind. 

Aims. The bulk of the Galactic CRs (GCRs) are accelerated by shock waves generated in supernova remnants (SNRs), a significant 

fraction of which occur in OB associations on a timescale of several 10^ years. We examine the effect of changing boundary conditions 

at the base of the galactic wind due to sequential SN explosions on the outflow. Thus pressure waves will steepen into shock waves 

leading to in situ post-acceleration of GCRs. 

Metliods. We performed hydrodynamical simulations of galactic winds in flux tube geometry appropriate for disk galaxies, describing 

the CR diffusive-advective transport in a hydrodynamical fashion (by taking appropriate moments of the Fokker-Planck equation) 

along with the energy exchange with self-generated MHD waves. 

Results. Our time-dependent CR hydrodynamic simulations confirm the existence of time asymptotic outflow solutions (for constant 

boundary conditions), which are in excellent the agreement with the steady state galactic wind solutions described by Breitschwerdt 

et al. (1991, A&A 245, 79). It is also found that high-energy particles escaping from the Galaxy and having a power-law distribution 

in energy (oc £"^') similar to the Milky Way with an upper energy cut-off at ~ 10'^^ eV are subjected to efficient and rapid post-SNR 

acceleration in the lower galactic halo up to energies of 10" - 10'^ eV by multiple shock waves propagating through the halo. The 

particles can gain energy within less than 3kpc from the galactic plane corresponding to flow times less than 5- 10* years. Since 

particles are advected downstream of the shocks, i.e. towards the galactic disk, they should be easily observable, and their flux should 

be fairly isotropic. 

Conclusions. The mechanism described here offers a natural and elegant solution to explain the power-law distribution of CRs 

between the "knee" and the "ankle". 

Key words. Galaxies: evolution - ISM: jets and outflows - Galaxies: starburst - supernova remnants - cosmic rays 



!^ 1 . Introduction power-law when going from massive to poor clusters of galax- 

, ies and groups. This has been interpreted in terms of a so-called 

It is now understood that galactic winds may form an im- "entropy floor" that dominates the gravitational heating of flie 
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; portant evolutionary stage in galaxy evolution. In particular in ^^^^^^j. ^^ j^e potential wells become shallower (PonmanHU 

the early universe, winds may be held responsible for pollut- ^^ jj^ ^ther words, the preheating of the IGM by starburst 

ing the intergalactic medium (IGM) with metals, possibly syn- driven galactic winds starts to dominate the release of gravi- 

thesized in Population III stars (e.g., Loewenstein, 2001| The ^^^^^^^^ energ y. In addition , the fairly high met aUicties of Z ~ 

discovei7 of Lyman-Break galaxies (e.g., .Steidel et al.t [19961 03 ^^ (^ g^ Mole ndi et all [19991: iKikuchiet ill [199^ found in 

[ Adelberger&Steideli [2000|) has given much support to this j^e intracluster gas may be explained by ejection of metals from 

idea. Evidence of winds at high redshift has been reported, galaxies by massive winds d uring starburst and also more quies- 

e.g., by Dawson et al. (2002), who serendipitously discovered a ^g„j phases (.Kanfereret all [2006.1 . 
strong asymmetric Ly„ emission line in a faint compact galaxy 

at z = 5.189, indicating an outflow velocity in the range of i^ the local universe, evidence for galactic winds mainly 

320 - 360 km s . stems from observations of starburst galaxies, such as M82 or 

In a recent study, it has been shown that many star- forming NGC 253. In these more evolved objects, starbursts are thought 

galax ies at 2 < z < 3 exhibit galactic outflows (Steidel et all to be triggered by a substantial disturbance of the gravitational 

l2010l) . Apart from metal enrichment of the IGM, galactic winds potential, e.g. by interaction with a companion galaxy. Finally, 

are also thought to preheat the gas. The X-ray luminosity ver- dwarf galaxies are prone to galactic winds because of their shal- 

sus temperature relationship (Lx - T) shows a deviation from a low gravitational potential well and their correspondingly low 

escape speeds, easily reached by supernova heated gas of a few 

Send offprint requests to: E. A. Doifi 10^ K. According to Bernoulli's equation, the typical outflow 
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velocity, if all thermal is converted into kinetic energy, is on the 
order of v ^ 250 km/s. In contrast, it is still widely believed that 
late-type spirals like the Milky Way should not exhibit winds, 
simply because the total thermal pressure in the disk is not suffi- 
cient to enable break-out through the extended gaseous HI layer 
and drive an outflow. 

The situation though may be different in very active star 
forming regions, where superbubbles locally inject a substantial 
amount of energy enabling hot gas to re ach considerable heights 
(lAvillez & Breitsch werd?. '2004', '2005"). A related feature was 
observed in Hq, by Reynolds et al. (2001) with WHAM. They 
find a giant loop extending out to ~ 1300 pc perpendicular to 
the galactic plane. The base of the loop is spa tially and kinemat- 
ically connected to the HI chimney found bv lNormandeau et alj 
(1199 6*). which is generated by the W4 star cluster As more en- 
ergy is deposited into this system by stellar winds and successive 
supernova explosions, we can expect the loop (the inner ionized 
part of the shell) to expand further and fragment due to Rayleigh- 
Taylor instabilities. The released hot plasma will then form the 
base of a galactic fountain, in which the gas may travel several 
kpc above the plane, but owing adiabatic and radiative cooling, 
thermal instabilities will progress and condensation of clouds 
(which ar e still grav itationally bo und to the galaxy) is supposed 
to occur (iBregmani [1980; Kahn, Il98lb lAvillezL l2000h . Only if 
the combined thermal energy of breaking-out superbubbles ex- 
ceeds the gravitational potential, a thermal galactic wind occur 
In normal spirals this will only happen locally, and more ex- 
tended winds are only expected in starburst regions like in M82 
orNGC253. 

However, as _ emphasize d by many authors in the past 



(e.g..lAxfordllT9"81.; Ipaviclil 119751; iBreitschwerdt et al.L I199U 
lEverett et al.L l2008l) . the above argument ignores the effect of 
cosmic rays (CRs), which have an energy density that is roughly 
in equipartition with the thermal gas and the magnetic energy 
densities. The fact that CRs escape from the Galaxy after a few 
tens of million years on average, setting up a pressure gradient 
pointing away from the disk, gives rise to the resonant excitation 
of MHD waves via the so-called CR streaming instability (e.g., 
[Kulsrud & Pearce, 1969). Thus CRs are coupled to the thermal 
plasma by scattering off the frozen-in waves, thereby helping to 
push the plasma against the gravitational pull. As a result, CR 
transport, which is mainly diffusive in the disk, will predomi- 
nantly become advective in the Galactic halo. Thus even in case 
of the Milky Way, a steady state wind flow could be set up with 
a global mass loss rate of ~ 0.3 Mpyr"^ and a low outflow base 
velocity of ~ 10 km/s ( Breitschwerdt et allll993l) . In fact, it has 
been shown that the diffuse broad-band soft X-ray background 
can be explained all the way from 0.3 keV to 1.5 keV by a su- 
perposition of emission from the Local Bubble and by a wind 
from the Galactic halo wit h both plasmas being in a s t ate of 
nonequilibrium i onization (IBreitschwerdt & Schmutzlerl 11994 
119991) . Recently, lEverett et al.l (|2008|) have shown fliat a kilo- 
parsec scale galactic wind of the Milky Way gives the best ex- 
planation for the longitude averaged soft X-ray emission near 
0.65 and 0.85 keV observed with ROSAT PSPC. Their model is 
further constrained by simultaneously fitting the observed radio 
synchrotron emission, resulting in a higher pressure and smaller 
disk surface area wind (lEverett et al.ll2010l) . 

Common to all investigations is the energy source due to 
more or less frequent supernova (SN) explosions, as well as to 
stellar winds, which contribute a minor fraction, and only during 
the initial phases, to the momentum and energy input, because 
only very massive O stars have sizable winds. For high stellar 
birth rates, therefore we can parameterize the results by the star 



formation rate (SFR) which relates the newly formed stars to the 
SN rate. In the case of low activity, the overall SN rate (types 
I and II) determines the properties of the global galactic wind. 
Assuming supernovae to be the sources of cosmic rays, we in- 
deed expect CRs to accompany outflows from the disk into the 
halo, unless their coupling is very weak, and diffusion is the 
only mode of transport. However, as CR observations suggest, 
the diffusio n halo only extends ~ 1 - 3 kpc perpendicular to 
the disk (see lBreitschwerdt et al.L 120021) . We therefore argue that 
CRs will always take part in the acceleration of an outflow, and 
although their effect may be small for starburst galaxies, it will 
dominate in the case of normal spirals. 

On the other hand, advection of CRs in the galactic wind 
flow may lead to reacceleration of escaping disk particles in the 
energy range 3 x 10'^ - 10'^ eV (i.e. between the "knee" and 
the "ankle") by propagating comp r ession al waves as has been 
suggested by IVolk & ZirakashviUl (|2004|) . These waves come 
from the interaction between fast outflows, originating in re- 
gions in the disk that have been compressed by the spiral density 
shock wave, hence have generated a faster CR driven outflow 
(cf. Breitschwerdt et al. 2002) and slower outflows from inter- 
arm regions. The steepening of the waves into moderate shocks 
of a sawtooth shape at large distances, where the flow is basically 
radial, leads to reaccelerating the particles close to the knee at a 
quasi-perpendicular shock (since the field is basically azimuthal 
there). Moreover, the supersonic galactic wind expands into the 
intergalactic medium with a given (but poorly known) pressure 
and must therefore be bounded by a termination shock, which 
itself can be the site of particle acceleration (Jokipii & Morfill, 
1 19851 Il987l) . It will, however, be difficult to observe these par- 
ticles since they are advected away from the Galaxy with little 
chance to diffuse upstream to the disk. Since MHD turbulence 
downstream of the ter mination shock is pre sumably very strong, 
it acts in the model of Volk & Zirakashvilii (2004) as a de facto 
reflecting boundary for the accelerated particles up to a certain 
energy. As the spiral density wave pattern is not in corotation 
with the disk and hence with the frozen-in halo field lines the 
associated forming halo shocks slip through the flux tubes. 

In this paper we present an alternative possibility of acceler- 
ating particles in a galactic wind flow beyond the knee. The basic 
idea is that individual star-forming regions, e.g. superbubbles or 
localized starbursts, have an energy output that is not constant 
in time. Moreover, their lifetime is typically less than 10** years 
and therefore smaller than the wind flow time. This will lead 
to time -dependent switch-off and switch-on effects of combined 
disk gas and CR outflows, which periodically push a piston into 
the gas. As we show in this paper, the galactic wind flow will be 
modulated by successive pairs of shocks (forward and reverse 
shock), where efficient particle acceleration can take place. Like 
the termination shock and unlike disk SNRs, these shocks will 
be long-lived, so the energy cut-off argument does not apply in 
this case because of finite size and acceleration time. They are 
thus strong candidates for generating the CR spectrum between 
lO''' eV and 10''^ eV by post-acceleration of GCRs, as we shall 
see. 

The paper is structured as follows. In Section |2] we describe 
the physical and numerical requirements of our model, and in 
Section|3]our general results from time-dependent galactic wind 
simulations for the Milky Way are presented, in particular the 
effect of particle diffusion and the propagation of shock waves 
within galactic winds. In Section|4|we introduce the interesting 
new possibility of accelerating cosmic rays by successive shock 
waves, which are driven by supernova explosions in the under- 
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lying disk, and propagate into the halo. A discussion and our 
conclusions are given in Section|5] 

2. Galactic wind model 

2.1. Physical equations 

When employing the usual physical notation we use the follow- 
ing time-dependent system of hydrodynamical equations supple- 
mented by three energy density equations for the thermal gas, 
£g, the cosmic rays, E^, and the Alfvenic waves £„: 



dp 

| + V.(pu) = 

^ + V-(pu:u) + V(Pg + P, + Pw) 
ot 

-HpVfSgal = 

dE„ 

— -^-hV-(£„u)-hP„V-u = r-A 

^ + V-((£e(u + rcVA)) + PcV-U 
Ot 



-VaVPc = V-(^V£e) 



at 



-hV-(£w(u + va)) + P„V-u 



h-vaVPc = 0. 



(1) 

(2) 
(3) 

(4) 
(5) 



The Alfven velocity is denoted by Va, the bulk gas velocity 
by H, its density and pressure by p and Pg, respectively, while 
the CR and MHD wave pressures are labeled by P^ and P„. 
The galactic gravita tional potential is Ogai, which we take from 
iBreitschwerdt et aP (Il991h . who used a Miyamoto-Nagai poten- 
tial for the disk and bulge, as well as a spherical dark matter 
component. Additional heating F and cooling A can be included 
and an energy sink for cosmic rays is incorporated by VaVPc, 
transferring energy from the cosmic rays to the waves due to 
the streaming instability. Energy sinks for the MHD waves other 
than adiabatic losses, e.g. ion-neutral or nonlinea r Landau damp - 
ing are not included (for a detailed discussion see lPtuskinlll997h . 
The former is neglected because the level of ionization can be 
easily maintained in a tenuous hot outflow plasma, and the lat- 
ter is expected to definitely become important, once the fluc- 
tuation amplitude reaches 6B/B ~ 1, when quasilinear theory 
breaks down. Since we start with a negligible wave field, assum- 
ing that all outgoing waves are generated by the resonant insta- 
bility, damping conditions in most model runs occur only farther 
out in the flow, where for example particle acceleration has al- 
ready taken place. However, in general the effect of nonlinear 
L andau damping shou ld not be neglected (see also discussion 
in lEverett et al.L 12008). Finally, in the case of thermally driven 
winds, such as in starburst galaxies, the contribution to the ther- 
mal pressure by wave damping is generally small, because the 
fraction of CR energy density is lower than the thermal one. 

The above system is closed by the equation of states relating 
the pressures to the energies by 



~ vyc,g,w ^ ) ^c 



(6) 



Assuming ultrarelativistic particles, we adopt y^ = 4/3, and tak- 
ing a monoatomic gas, we set y^ = 5/3. A value of yw = 3/2 
corresponds to a wave energy density E„ - {{6B)^)/4-n of the 
mean squared fluctuating magnetic field component. Since the 
magnetic wave pressure is given by P^, = (((5B)^)/87r, we ob- 
taina value of -y^ = 3/2 by analogy. We also have t o spec- 
ify the mean diffusion coefficient k for the cosmic rays (iDrurvi 



Il983h . which is averaged over the particle distribution func- 
tion (in space and momentum), and is usually taken to be a 
free parameter. According to cosmic ray propagation models 
fe.g.. iGinzburg & PtuskinL [l985l) . k is estimated as an order of 
10^9 cm v. 

In the self-consistent CR propagation model of Ptuskin et 
al. (1997), the field parallel diffusion coefficient is calculated 
from the MHD turbulence level arising from CR streaming, 
with generation balanced by nonlinear Landau damping, so that 

^11 - 10^^ (^7) cm^s"', for Galactic halo conditions, where 
^ ^ 4 is the power-law index in the particle momentum distri- 
bution near CR sources, and p,m, and Z are the particle mo- 
mentum, mass and nuclear charge, respectively. The influence 
of particle diffusion and acceleration is discussed in Sect. 13.11 
and mor e detailed properties of the system of Eqs. ([T}-(|6) can be 
found in IBreitschwerdt et all (1 199 ih . 

These equations are solved in a ID flux tube geometry, where 
all quantities depend on the projected distance z from the galac- 
tic plane and the time t. To fix ideas we specify the geometric 
properties of the flux tube, which represents a transition from 
planar to spherical flow, and is in its simplest form written as 



A(z) = Ao 



ZO 



(7) 



where zo defines a typical scale above which the flux tube opens 
due to spherical divergence (typically of the order of a galac- 
tic radius). For short distances z <K zo, the flow is essentially 
planar and as seen e.g. in Table [T] the lower values of zq lead 
to faster shocks propagating along the flux tube. Aq denotes the 
cross section at our reference level of z = 0, which cancels out 
in all equations. 

It has been emphasized by some authors that the equa- 
tions should be w ritten and solved in general ellipsoidal 
( Fichtner et al.Lll99L) or at least in oblate spheroidal coordinates 
(IK0III991]) . A test of the simple formulation of Eq. (O in com- 
parison to oblate spheroidal coordinates has shown that the max- 
imum deviation is an enhancement of the mass loss rate by about 
15% for a flux tube at a Galactocentric distance of Ro - 10 
kpc (iBreitschwerdtl Il994l) . The discrepancy even decreases for 
shorter distances, because the effect of a radial (as opposed to 
vertical) flow component becomes progressively weaker In ad- 
dition, it is more likely that the injection of supernova heated 
gas into the base of the halo also imparts a z-component of mo- 
mentum to the flow, again enhancing the vertical component. 
Thus, for practical purposes, Eq. (|7} seems fairly robust and 
time saving, and only a complete solution of the axisymmetric 
MHD equations, in which the surfaces of magnetic pressure are 
calculated self-consistently by a Grad-Shafranov type equation, 
marks a substantial improvement. Here we are primarily con- 
cerned with the interplay of various physical processes and the 
role of cosmic rays as a driving agent, and so we take advantage 
of the simple but flexible flux tube formulation in Eq. Q. 

As seen e.g. in Eq. ( fTSl l, the assumed flux tube geometry also 
determines the radial structure of the magnetic field. However, 
edge-on observations of spiral galaxies like NGC 5775 reveal 
large sc ale magnetic fields of similar topology in their galac- 
tic halo (iSoida et al.L l201lh . Also rec ent 3D-numerical simula- 
tions of a cosmic ray driven dynamo (iKulpa-Dvbel et all 1201 ll) 
show that random magnetic fields will evolve into large-scale, 
ordered magnetic structures that reveal a quadrupole-like con- 
figuration with respect to the galactic plane. Observationally, as 
well as theoretically, assuming a prescribed flux-tube geometry 
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thus seems appropriate for simplifying the time-dependent evo- 
lution of galactic winds through this ID-geometry. 

To summarize the previous section we have shown how to set 
up a system of equations that describes a galactic outflow driven 
by the interplay of the gravitational potential and the interaction 
of thermal gas, cosmic rays, and waves in an ad-hoc prescribed 
geometry. Several outflow simulations are discussed in detail in 
the following sections. The time-dependent effects of dissipa- 
tion, heating and cooling, as well as particle acceleration lead to 
a number of features, which in most cases can be followed until 
a time-independent, stationary solution is reached. 

2.2. Boundary conditions 

Boundary conditions close to the disk have to be specified, al- 
lowing for a transition from a subsonic to a supersonic out- 
flow. In steady state this is only poss ible if the outw ard pressure 
tends to zero at infinity, as shown bv lParkerl(ll958b for the solar 
wind. The terminal wind velocity can be easily estimated from 
Bernoulli's equation. At the inner boundary located at a refer- 
ence level, we have to fix, e.g., the gas density, the gas, and wave 
pressure. The cosmic rays are determined by an outward directed 
flux 



2.3. Numerical metliod 



F(. n — 



Tc 



(uo + Vao)^c,0 - 



dP, 



Tc - 1 dz 



(8) 



where the first part describes the advected particles through the 
combined wave speed uo + Vao, and the second term represents 
the diffusive flux across the inner boundary. Since energetic par- 
ticles always leave the galaxy within a so-called residence time 
(about ~ 10^ yrs for the Milky Way), we have to ensure Fc,o > 0. 
If we use a typical length scale L to approximate the diffusive 
term in Eq. ^ by Pcfi/L, we can calculate on which scale the 
cosmic ray flux would vanish, i.e., F^fi = OorL ~ /^/^^(uo + Va). 
Even in the case of low k - 10^^ cm^s"', we obtain L > I Mpc, 
hence F^fi > for all realistic galactic parameters since the 
Alfven velocity is |vaoI = 6.9km/s for our initial conditions. 
The last quantity to be specified at the inner boundary is the ve- 
locity Uo, which in case of a steady-state solution is determined 
by the conditions at the critical point, thereby fixing the mass 
loss rate. 

At the outer boundary we impose simple outflow conditions; 
i.e., no gradients should exist at large distances from the galactic 
plane, typically at z = lOOOkpc. As studied in various compu- 
tations the details of the outer boundary do not change the over- 
all behavior of the solutions. In steady-state flows, the domain 
of influence of the inner boundary conditions only encloses the 
subsonic flow region. The location of the critical point, how- 
ever, is not completely independent of the outer boundary con- 
ditions. In particular, it is required that the pressure at the sonic 
point exceeds the pressure at the outer boundary. Otherwise, no 
transsonic flows are possible and only breeze solutions will ex- 
ist. For time-dependent flows such restrictions do not apply in 
general, and only in the time-asymptotic limit do comparisons 
with steady-state flows make sense. 

As described in the following sections, we have computed 
time-dependent galactic winds arising from changes of the inner 
boundary conditions. In particular, we studied how these varia- 
tions can be traced in the outflow and how the perturbed solution 
evolves towards a stationary flow. 



Since the numerical method has been described in detail (iDorfil 
1998'), we emphasize here only some important features. First, 
all equations are cast into a volume-integrated form to ensure 
conservation of global quantities such as mass, momentum, and 
energy. Such a finite volume formulation is also suited to an 
adaptive grid that concentrates grid points to resolve steep gra- 
dients. Second, the equations are discretized in an implicit way 
so that for stationary solutions the time step exceeds the typi- 
cal flow time by several orders of magnitude. Third, an artificial 
viscosity with a variable lengt h scale in the tensor formulation 
(iTscharnuter & Winkleri 1 19791) is adopted to broaden the shock 
waves over a finite number of grid points. However, since we 
also use an adaptive grid, the typical broadening length scale 
/vise can be specified to be smaller than all physical length scales, 
which enables correct computation of particle acceleration; i.e., 
'vise ^ ^/ms with Ms denoting the shock speed. The first-order 
Fermi mechanism requires that the energetic particles are ex- 
posed to a sharp discontinuity to gain energy by being reflected 
between the upstream and downstream media. 

Applying the adaptive grid to these computations has the ad- 
vantage that the position of discontinuities can be located very 
precisely through a clustering of grid points without specifying 
the type of nonlinear wave beforehand. The propagation of such 
nonlinear waves is plotted, e.g., in Fig.|4] which enables a direct 
comparison with analytic estimates. 

Computations based on an implicit method have the advan- 
tage of allowing calculations of stationary solutions with the 
same numerical method as for the time-dependent solutions. 
As a result several solutions with fixed inner boundary val- 
ues develop towards the stationary solutions, all of which can 
be reached from different evolutionary paths. We found from 
such a large number of calculations that the stationary solutions 
look rather stable with respect to perturbations and that time- 
dependent flows relax to the steady-state after the initial pertur- 
bations have propagated through the computational domain. For 
our Milky Way with a typical flow speed around 500 km/s we 
get a typical relaxation time of about lO** years for a spatial re- 
gion of 1 Mpc. The initial mo dels can be taken directly fr om the 
time-independent solutions of lBreitschwerdt et al.l(ll991l) . 

The actual discrete equations are then solved by a nonlin- 
ear Newton-Raphson procedure and are given in the appendices. 
The total number of spatial grid points is set to 300 throughout 
the computations. Due to our six physical variables (gas den- 
sity, gas velocity, thermal pressure, cosmic ray pressure, wave 
pressure and radial position) at every grid point, each time step 
requires 1800 new unknowns to be computed. 

3. Results 

3. 1 . Effects of cosmic ray diffusion 

Before we compute time-dependent solutions of the system of 
CR-hydrodynamics in a flux tube geometry (Eqs.[T]|5) we inves- 
tigate the effects of cosmic ray diffusion in steady state solutions. 
The additional diffusive term written as 



V-(^VPc) 



(9) 



increases the order of the system of equations and no simple and 
straightforward analysis of critical points can be used to charac- 
terize the solutions. As emphasized earlier the implicit character 
of our numerical methods makes it easy to implement such an 
additional transport mechanism and in Fig. [T]the influence of a 
mean diffusion coefficient k is shown. 
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Fig. 1. The time-asymptotic flow structure of a galactic wind up to the innermost 200 kpc with different cosmic ray diffusion 
coefficients Tc. The case k - lO^^cm^s"' (fufl line) is almost identical to /<• = 0, and noticeable differences occur only for Ic - 
10^'' cm^ s"' (dotted line), k - 3-10^^ cm^ s ' (dashed line), and k - 10^*' cm^ s ' (dashed-dotted line). The different panels represent: 
a) gas density, b) gas velocity, c) thermal gas pressure, d) ratio of gas pressure to cosmic ray pressure. 



The spatial wind profiles plotted in Fig.[T]exhibit the changes 
of the overall flow structure caused by different cosmic ray dif- 
fusion coefficients in the astrophysical relevant parameter range 
qIk- 10^^ - lO-'" cm^s"' within the innermost 200 kpc. Starting 
at a reference level of 1 kpc, we follow the flow up to a distance 
of 1000 kpc in order to reach asymptotic values. For these com- 
putations the numerical parameters correspond to our Galaxy, 
and we keep fixed most quantities at the inner boundary, i.e. the 
gas density of po = 1.67- lO^^^gcm"-', the thermal gas pres- 
sure fg_o - 2.8 ■ 10^'^dynecm"^, the vertical magnetic field 
component Bq = 1 y^G, and the initial wave pressure Pw.o - 
4- 10^'^dynecm"^, which corresponds to magnetic fluctuations 
at a low level of 0.1 Bq. In the case of no diffusion (k - 0), the 
mass flux and therefore the initial velocity mq are determined by 
the critical point (CP) conditions and left as a free inflow bound- 
ary for the diffusive cases. To investigate the effect of cosmic 
ray diffusion we fixed the total cosmic ray flux Fc,o at the inner 
boundary through condition (jSJ. 

Since the diffusion length scale L associated with a flow 
velocity U and diffusion coefficient k is determined through 
L ^ kIU , i.e., for a typical value of L'^ ^ 300km/s and of 
K ^ lO^^cm^s"', we get L ^ 1 kpc, so significant modifications 
of the large scale behavior are only expected for values of the 
mean diffusion coefficient k > lO^^cm^s"', as can be seen in 
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Fig. 2. Relative mass-loss rates (filled symbols) M in units of 
Mo, the mass-loss rate at ^ = 0, and corresponding terminal ve- 
locities (open symbols) in [km/s] for a flux tube located in our 
Galaxy at a galactocentric distance of 7?o - 10 kpc as a func- 
tion of the mean cosmic ray diffusion coefficient k. Only values 
of /? > lO^^cm^s"' reduce the mass loss rate significantly, but 
also increase the final wind velocities owing to lower thermal 
gas ejections (see text for details). 
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Fig.[T]for the Milky Way. Since the energetic particles can make 
an important contribution to the overall galactic gas outflow, the 
coupling of CRs to the gas is an important process that alters the 
mass loss rate of the wind. In the case of diffusive losses, CRs 
can escape from a galaxy without having to drag thermal gas 
along. To compare the solutions, the total cosmic ray flux Fc,o 
is kept constant in Fig. [T] and therefore an increase in the diffii- 
sive losses has to be compensated by a corresponding decrease 
in the advection of particles at the reference level; i.e., uq de- 
creases. Since the mass flux (per unit area) is given by pquq, dif- 
fusion reduces the amount of thermal material within a flux tube. 
Consequently, any increase in the diffusion coefficient k leads to 
lighter and faster winds with lower mass loss rates. Since the 
parameter s for this computation are c hosen for our own galaxy 
similar to iBreitschwerdt et all (Il991b . cosmic rays provide the 
dominant pressure contribution above lOkpc for all values of k 
as can be inferred from Fig. [TJl. This is because even for the 
highest value of k the diffusion length does not exceed lOkpc. 
We emphasize that the solutions shown in Fig. [T] for k = for 
gas density, outflow velocity, thermal, CR, and wave pressures, 
all as a function of distance, are i n excellent agreement with the 
steady-state solutions obtained in IBreitschwerdt et al.l(ll99lh . 

The cosmic rays affects the inner boundary through the in- 
going cosmic ray flux F^fi (Eq. [8]) as well as the transport of 
material through the flux tube thanks to the available cosmic ray 
pressure gradient. The increase in the mean cosmic ray diffusion 
coefficient k leads to changes in the mass-loss rate and the final 
velocity as is shown in Fig. [2] Since in the adopted flux tube 
geometry for stationary winds the mass-loss rate is given by 



M - Aop()M() - nR pquq , 



(10) 



we have plotted in Fig.|2]the mass loss relative to the case with- 
out cosmic ray diffusion, i.e. k - Q, and denote this value by 
Mq. Therefore we get rid of the scaling by the initial cross sec- 
tion of the flux tube Aq = nR^. This figure shows how the total 
mass loss is decreasing as the relative contribution of cosmic 
ray diffusion oc kdPc/dx increases compared to the advected 
flux oc (mo + va)Pc,o for a fixed value of the total CR-flux F^s) 
at the inner boundary (cf. Eq. |8]l. As stated before, the typical 
length scale of cosmic ray interactions on a flow with velocity 
U is given by k/U, which alters the cosmic ray pressure gradi- 
ents available to drive the outflow. As seen in Fig. |2] values of 
k i lO^'^cm^s"' only have a small effect on the outflow char- 
acteristics, because such values are associated with length scales 
around 1 kpc taking a typical flow velocity of U ~ 300km/s. 
Comparing the results for various simulations of galactic winds 
from our Galaxy (Fig. [T}, we conclude that only mean cosmic 
ray diffusion coefficients above k k. 10^'' cm^ s"' affects the final 
outflow velocities and mass loss rates significantly. Constraints 
based on galactic cosmic ray observations indicate a diffusion 
coefficient on the order of ~ 10^'^ cm^ s ', for nucleons below 
about a GeV (e.g., Axford, 198 1 ), which represent the bulk of the 
particles, and somewhat higher values ~ 10""' cm^ s"' for parti- 
cles with higher energies. In addition, the generation and prop- 
agation direction o f Alfven waves are influenc ed by the cosmic 
ray gradients ("e.g.. 'Breitsch werdt et al.l Il99ll) . However, if pa- 
rameters that are appropriate for different galaxies only make a 
minor contribution to the wave pressure compared to the gas and 
cosmic ray pressures, this effect is almost negligible here (see 
also Figs. Hi and|5jl). 

All results reported here were calculated for a Miyamoto- 
Nagai galactic potential as used for our Milky Way at a galac- 
tocentric distance of Rq = 10 kpc. This is slightly higher than 
the lAU value for the solar circle, but has been adopted here 



in order to compare results with IBreitschwerdt et aP (Il99lb . 
The dependence on the location within a galaxy enters only di- 
rectly through a variation in the galactic potential Ogai, and in- 
directly by changes i n gas density, pressure etc., as discussed in 
IBreitschwerdt et al] (|I99I|) . Integrating their galactic winds so- 
lutions over the whole galactic disc yields a total mass loss rate 
of Mgai ~ 0.3 Mq yr"'. This can serve as an upper limit of the 
global galactic mass loss rate since the effect of cosmic ray dif- 
fusion (Fig.|2| will reduce this value, which was calculated from 
purely advective solutions. 

Since the outflow velocity in case of thermally driven winds 
is usually more than an order of magnitude above the escape ve- 
locity, the mode of CR transport is immaterial here, since the 
mass loss rate itself is mainly determined by the thermal pres- 
sure. This is a consequence of the softer equation of state of the 
CRs, as compared to the thermal plasma, which leads to an effi- 
cient lifting of the halo weight at larger heights, when the mass 
density is akeady smaller In terms of galactic wind flows, the 
thermal plasma gradient lifts up the gas close to the disk, and 
the CRs act as a kind of afterburner, accelerating the sluggish 
thermal winds to higher velocities. This effect is obviously exac- 
erbated by an increasing CR diffusion coefficient. 

The flow time through the wind is given in our computational 
domain by 



fflo 






dz 
<i) 



(11) 



where the boundaries have to be taken e.g. between 1 kpc and 
1 Mpc. In accordance with higher outflow velocities for increas- 
ing cosmic ray diffusive transport the flow times up to 1 Mpc 
decrease from 4.1 ■ 10^ years (k - 0) to 2.8 ■ 10^ years in the 
case of k - 10^"cm^s"'. The flow times or final velocities do 
not scale directly with the mass loss rates (Fig.|2]i since increas- 
ing the diffusion k changes the inner boundary conditions (Eq. |8]l 
and also lowers the cosmic ray gradients needed to further accel- 
erate the outflow. The sonic point of the galactic winds shrinks 
from 108 kpc to 48 kpc. 

3.2. Time-dependent flow features 

Assuming that the flux tube geometry is given and fixed in time 
we expect at least time-dependent inner boundary conditions. In 
a realistic case of repeated SN-explosions or a starburst event (on 
larger scales), the inner variations will be determined by the tem- 
poral history of such explosive events. Every change at the base 
of the flow introduces additional disturbances, which then prop- 
agate outwards and modify the flow conditions along the flux 
tube. To keep these structural changes simple and to illustrate 
the time-dependent behavior of a galactic wind, we have chosen 
a stationary model as our initial model with a flux tube opening 
distance zo = 15 kpc and a mean cosmic ray diffusion coefficient 
ofk- 10^^ cm^/s, plotted also as dotted line in Fig.[T] The prop- 
agation of waves is driven by a sudden increase in gas and cos- 
mic ray pressure at the inner boundary located at 1 kpc, where we 
have raised the pressures by a factor of 10 over a time interval of 
10*' years. All other parameters are kept constant. Since we use 
an implicit numerical method, we can follow the propagation 
of these pressure waves through the whole computational do- 
main until a new stationary solution has been established within 
1 < z/[kpc] < 1000. Typicafly, the flow time (Eq. [TTJ up to 
1 Mpc is of a few 10^ years. 

As plotted in Fig. |3]for the innermost 100 kpc the changes 
at the inner boundary result in a propagation of strong nonlinear 
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Fig.3. The time-dependent flow structure of a galactic wind up to 100 kpc at six different times at 3.4-10^, 2.0-10^, 3.6-10^, 5.4-10^, 
7.7-10^ and 1 .0-10^ years. At f = the gas and cosmic ray pressure has been increased by a factor of 10 at the inner boundary of the 
flux tube, simulating a violent SN-explosion. Panel a) plots the development of the density structure where the forward as well as 
the reverse shock become visible. The velocity in [km/s] is depicted in panel b) where the particle acceleration broadens the shock 
fronts in time. In panels c) and d) the curves correspond to the gas and the cosmic ray pressure (solid line) together with the wave 
pressure (dashed lines), respectively. The velocities of these features are also plotted in Figs. |4]and|6] (see text for more details). 



waves along the flux tube since the conditions at the base of a 
galactic wind have been altered dramatically. In the simulations 
reported here, we keep the situation as simple as possible and see 
how pressure waves are initiated and gen erate a flow p attern sim- 
ilar to a supernova remnant (SNR) (e.g.. lDorfilll991 ). A strong 
forward shock wave propagates through the unperturbed wind 
structure, accompanied by a reverse or terminal shock, which 
slows down the newly developed outflow. However, owing to 
thehigh density and pressure gradients within the initial flux 
tube, the overall properties of the shock waves cannot be charac- 
terized by a rather simple self-similar behavior, which is typical 
of the SNR case. At both shock waves, particles are accelerated, 
which modifies the shock structure itself, and the dissipation of 
waves heats the thermal gas. The maximum momentum p^ax of 
the accelerated particles attainable in such a shock configuration 
is plotted in Fig. |7] resulting in a reacceleration of the Galactic 
cosmic rays already driving the outflow. At the moment we have 
not included any radiative losses of the thermal plasma in our 
computations. Radiative cooling should not be dramatic close to 
the disk where the temperature is still above 10^ K and where the 
bulk of the particles are post-accelerated. At greater distances. 



when the temperature and density have dropped, cooling times 
(which are proportional to 1 /p) still remain high. 

In Fig. [3^ the gas density clearly exhibits the separation of 
the flow structures between 3.6-10^ years and 1.0-10** years since 
the forward shock propagates finally at a speed of 930 km/s, 
whereas the reverse shock stays behind at 840 km/s. Already 
after 10' years the velocity difference has accumulated to a spa- 
tial separation of about 1 kpc. The variation in these velocities 
caused by diff'erent flux tube opening scales zo (see Eq. |7]i is 
given in Table[T] The gas compression ratio at the forward shock 
varies in time (cf. Fig. HJi since the flow is strongly affected by 
the cosmic ray pressure contributions having an adiabatic index 
of Tc = 4/3 as well as by the acceleration of cosmic rays, which 
tends to broaden the shock transition (cf. Fig. |4]for more details). 

The gas velocity portrayed in panel b) of Fig. |3] exhibits the 
importance of the forward shock where the outflow speed within 
the innermost 100 kpc increases from values less than 200 km/s 
to more than 800 km/s. Although visible in the gas density and 
gas pressure, the reverse shock plays a minor role for the thermal 
gas during this initial expansion since the flows quickly acceler- 
ates down the large gradients. Only at distances z ^ 200 kpc does 
the Mach number of the reverse shock increase again when the 
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Fig. 4. The locations of the two shock fronts (forward shock: 
filled squares; reverse shock: open diamonds) for galactic dis- 
tances z < 100 kpc and compression ratios of the gas as a func- 
tion of time. Already at early times (f < 20- 10''years) the shock 
propagation almost becomes constant. The solid line denotes 
the compression ratio of the forward gas subshock, the dashed 
line depicts the total compression ratio including the contribu- 
tion form the cosmic ray precursor The dashed-dotted line plots 
the gas compression ratio of the reverse shock growing in time 
(see text for more details). 



flow propagates in background with small gradients. The spa- 
tial variations in the post-shock region are caused by gas pres- 
sure gradients (s. Fig. |3];), which are not smoothed by diffusive 
transport processes, such as the corresponding cosmic rays in 
Fig-EJl, but are affected by additional heating from the dissipa- 
tion of waves. Due to the pressure increase at the base of the 
flux tube the gas velocity at the inner boundary changes from 
uq - II km/s to a new asymptotic value of 5 17 km/s. Since the 
density po is kept constant, the mass loss M oc po^o also in- 
creases by a factor 517/11 ~ 50 from this flux tube into the 
intergalactic space. Clearly, all the above values depend on the 
galactic potential, as well as on the flux tube opening zq (see 
Table [U. As a result, these simulations may serve as a typical 
example for the case of our Galaxy. 

The features in the gas pressure are dominated by the for- 
ward shock as shown in Fig. [3}; within the innermost 100 kpc. 
Since the flow is accelerating down the strong density gradient, 
the Mach number of the forward shock increases. As mentioned 
before, the reverse shock plays a minute role for the thermal gas 
at these stages as seen, e.g., in the rightmost curve by a small 
pressure jump at f = 4.8 ■ 10^ years, at a distance of 62 kpc. 

The cosmic ray pressure (Fig. |3]l) dominates the flow struc- 
ture as can be easily inferred by comparison with the gas pres- 
sure of panel c). Energetic particles contributing to the cosmic 
ray pressure are mainly accelerated at the forward shock, and 
the shock structure is smoothed when the precursor region de- 
velops. Since we have adopted a mean cosmic ray diffusion co- 
efficient of K = 10^'cm^/s the typical acceleration time scale 
yields t ~ k/u^ ~ 1.5- 10^ years for an initial shock velocity 
of 450 km/s. The spatial scale corresponding to this accelera- 
tion time is given by / ~ k/us ~ 0.8 kpc, which is very close 
to the inner boundary. As also discussed in Sect. |4]and depicted 
in Fig. |7]the particles are accelerated rapidly and close to the 
galactic plane. Although the reverse shock does not play an im- 
portant role for the overall cosmic ray pressure at the early out- 



flow stages (see also Fig. |4]l a small number of individual par- 
ticles will gain additional energy from this second shock wave 
(see also Fig|7]i. 

In Fig. [3}J we have also plotted the radial variations in the 
wave pressure P„ by dashed lines, where the corresponding 
wave energy equation (|5j contains the dissipation term vaVP^ 
heating the thermal gas. The typical shape of these curves clearly 
exhibits the location of the forward and the reverse shocks, seen 
as pronounced changes of P„. We can also infer the impor- 
tance of particle acceleration by the smoothed edges. Again, 
the wind parameters typical of our Milky way show the impor- 
tance of the excited MHD-fluctuations P„ - {{6B)^}/?,n with 
^w - (Tw - l)£^w and yw = 3/2 for driving the outflow. Without 
the additional diffusive transport due to V ■ (kVE^) acting on the 
cosmic ray pressure P^ the wave pressure P„ is largely affected 
by the reverse shock seen, e.g., at z = 73 kpc for the rightmost 
plotted time of f - 1.0- 10^ years. Within the zone of the two 
shock waves the wave pressure P^ remains above the thermal 
pressure Pg. However, we do not want to stretch this argument 
any further, because some care must be applied to the impor- 
tance of wave pressure. As mentioned before, nonlinear Landau 
damping will limit wave growth severely, resulting in additional 
thermal heating of the plasma. 

Since the development of shock waves plays an important 
role in re-accelerating the existing population of Galactic cosmic 
rays (see Sect.|4]i, we explored the propagation of these flow fea- 
tures in Fig. |4] within the first 10** years after our initial pressure 
increase at the inner boundary simulating multiple SN explo- 
sions. After the initial change, the inner boundary conditions are 
kept constant, causing nonlinear waves moving outwards at dif- 
ferent speeds, leading to a growing separation between the flow 
features as can be seen in Fig. |4] The velocities can be directly 
estimated from the slopes and are given in the reference frame 
of the galaxy and reach 930 km/s for the forward shock, and 
840 km/s for the reverse shock, respectively. Unlike the solar or 
the galactic wind, where the termination shock is almost station- 
ary with respect to the source frame, or in a supernova remnant, 
where the reverse shock is moving inwards thereby heating the 
ejecta, the reverse shock here is largely convected outwards by 
the galactic wind flow. In other words, its relative velocity with 
respect to the contact discontinuity (~ 80 km/s) is much too low 
compared to the galactic wind speed of ~ 900 km/s, to propa- 
gate back to the disk. This is a direct consequence of the density 
profile in the undisturbed galactic wind medium, falling off as 
1 /z^ as we shall see. The shock velocities are constant, which 
can be easily understood in the framework of similarity solu- 
tions. If we disregard any processes containing explicit length 
and time scales during the formation of the flow, the propagation 
of a fluid element in an ambient medium with density decreas- 
ing by a power-law, p oc z^, is similar to a stellar wind type 

flow, in which the outer shock travels a distance Rs oc t^*i' . Since 
for an SN or starburst driven outflow, the terminal wind speed 
is reached very quickly, and therefore the steady-state wind den- 
sity behaves as p„, oc (m^A(z))"', where A(z) is given by Eq.|7] At 
some distance from the plane, A(z) oc z^, therefore /3 = -2, and 
thus /?s oc f or the shock velocity m, = dRJdt ~ const., clearly 
visible after f ^ 2- 10^ years or z ^ 12 kpc for our Milky Way pa- 
rameter in Fig.|4](filled squares). The same must be true for the 
reverse shock (open squares in Fig. 0, which sees a constantly 
decreasing density profile, again oc l/z^, as it tries to propagate 
inwards, but is effectively convected outwards. This is also, by 
the way, the reason the reverse shock remains strong for a long 
time (~ 10^)years, because it does not run into medium with 
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Fig. 5. The time-dependent flow structure of a galactic wind up to 1 Mpc at seven diff'erent times at 3.7-10^, 1.1- 10*^, 2.3-10**, 3.7-10**, 
5.3 - 10**, 7.7- 10^, and 1.0- 10^ years. At f = the gas and cosmic ray pressure has been increased by a factor of 10 at the inner 
boundary of the flux tube, simulating a series of violent SN-explosions within a short time interval. In Panel a) the development of 
the density structure is plotted, where both the forward and the reverse shocks become visible. The velocity in [km/s] is depicted in 
panel b) where the particle acceleration broadens the shock fronts in time. In panels c) and d) the curves correspond to the gas and 
the cosmic ray pressure, together with the wave pressure (dashed lines), respectively. The velocities of the propagating features are 
also plotted in Figs.|4]and|6](see text for more details). 



increasing density, where it would have to compress and set an 
increasing amount of mass into motion. However, the diff'usive 
transport of cosmic rays, as well as the dissipative heating of 
the thermal gas is important during the early expansion phases 
and the strength of the reverse shock is therefore best visible in 
the wave pressure (Fig. |3]J, dashed lines). This means that the 
reverse shock also remains a site of efficient particle accelera- 
tion, although particle escape would be primarily away from the 
galaxy similar to the galactic wind termination shock. 

The solid line in Fig. |4] represents the compression ratio of 
the gas subshock within the forward shock (filled squares), also 
visible in Fig. [3^. The shock structure also exhibits a cosmic 
ray precursor where the cosmic rays diffuse ahead of the wave 
thereby reducing the incoming plasma speed and converting adi- 
abatically the kinetic energy into particle pressure. The thermal 
gas is compressed in this precursor region, and the overall com- 
pression ratio, including this effect, is plotted in the dashed curve 
of Fig.|4] leading to compression ratios in excess of 4. The varia- 
tion in the compression ratios is due to the changing background 
of the initial flow structure (cf. Fig.[T] dotted lines), which is ba- 
sically confined to the innermost 50 kpc. The dashed-dotted line 



depicts the gas compression ratio of the reverse shock, which 
increases in time, since the wind velocity at the inner bound- 
ary rises, and the flow has to be decelerated to the almost con- 
stant post shock velocity of the forward shock. We emphasize 
that both shock waves travel at different speeds, as summarized 
in Table [T] and that the most energetic particles can pick up the 
various velocity differences to get reaccelerated (see also Fig.|7). 

In Fig. |5]the variables are plotted between the inner bound- 
ary at 1 kpc and the outer boundary located at 1000 kpc at 
seven different times. The first time step depicted corresponds to 
3.610^ years, followed by 1.1-10**, 2.3-1 0^ 3.7-1 0^ 5.3-1 0^ 7.7-1 0^ 
and 1.03-10^ years. The density, velocity, and pressure structures 
in Fig.|5]clearly exhibit both shock waves and the contact discon- 
tinuity, which is more difficult to localize, since several pressure 
contributions are important within this zone. The slowdown of 
the gas velocity of about 100 km/s is caused by the smoothed 
reverse shock (Fig.|5]3). Both shock waves energize cosmic rays 
and in particular the further acceleration of particles at the re- 
verse shock is clearly visible through the local maxima in the 
corresponding P^-curves of Fig.|5}l- The increase in the particle 
pressure, as well as the wave pressure (dashed lines) reveals that 
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Fig. 6. The various velocities occurring in the time-dependent 
solutions also plotted in Fig. |5]as a function of time. All veloc- 
ities are given in the rest frame of the galaxy. The solid line de- 
picts the velocity of the forward shock v, and the dashed line the 
reverse shock v^. The highest gas velocity (filled squares) gives 
the upstream velocity of the reverse shock V2- The open triangles 
plot the post shock velocity vi of the reverse shock. The open 
diamonds correspond to the gas velocity vo three diffusive scales 
k/us in front of the forward shock. In less than about 7-10^ years 
or r < 50kpc, all the velocities become constant and are also 
given in Table [T] 



the reverse shock remains a site of efficient particle acceleration. 
Since after t ^ 1.4- 10^ years, the reverse shock has reached a 
distance of more than 100 kpc, the particles would tend to move 
away from the galaxy as at the galactic wind termination shock. 
In the forward shock the outflowing gas is pushed from 
400km/s to 800km/s. Consequently, the overall flow time 
through the wind up to 1 Mpc is shortened from the initial model 
of 3.2- 10^ years to 1 .08-10^ years corresponding to a factor of 3. 
The latter timespan is also identical to the time the initial pertur- 
bation needs to travel to the outer boundary located at 1 Mpc. 

4. Particle acceleration in galactic winds 

As seen in the previous section, a temporal variation due to 
changes in the star formation rate (and hence SN-rate) leads to 
an increase in the gas pressure at the bottom of a flux tube. Such 
time-dependent boundary conditions at the reference level re- 
sult almost inevitably in (several) shock waves propagating into 
the halo (and possibly coalescing) as outlined in the previous 
section. Under such conditions (two shock waves separated by 
a contact discontinuity, cf. Fig. |5]l energetic particles can gain 
energy through the first-order Fermi-mechanism. The bulk of 
the particles, which enter the diffusive shock acceleration pro- 
cess, are already energized CRs from the disk sources so that 
there is no need to worry about injection. These particles in- 
stead undergo a reacceleration process. In contrast to the com- 
pression waves of the "slipping interaction regions" discussed 
inlVolk & Zirakashvili (2004), the shock waves considered here 
can be strong. Since both forward and reverse shocks are run- 
ning down a density gradient, they can strengthen considerably 
already close to the disk (see Fig.|5]l. 

Diffusive particle transport along the magnetic field lines 
of the flux tube in the time-dependent galactic wind flow thus 
opens the possibility to accelerate cosmic rays to energies well 
beyond the knee, if large-scale shock waves propagate through 



the wind. Ad opting for simplicity the test particle picture (e.g. 
iDrurvi Il983h we can integrate the equation 



dp„ 



dt 



(12) 



to estimate the maximum momentum /?max reached during the 
life time of a shock wave. We use the particle momentum in units 
of inc. The diff'usive acceleration time scale at a shock wave with 
upstream velocity mi and downstream velocity U2 for a planar 
shock (the radius of curvature is much larger than for supernova 
remnants) is given by 



i Kl K2 

race = ~ + ~ 

M] — U2 \ Ml M2 



(13) 



The diffusion coefficients a-i 2 - k(p) depend on the particle 
momentum p. If fully developed turbulence is assumed the so- 
called Bohm limit, where the mean free path is essentially re- 
duced to a particle gyroradius, can be adopted, and for a particle 
with mass m and charge Ze the diffusion coefficient reads as 



1 mc^ 



3ZeB .JTTpi 



(14) 



The Bohm limit means that the field is completely random on a 
scale of the order of the gyroradius, and thus represents a min- 
imum diffusion coefficient. It is indeed questionable whether 
higher energy particles find a wave field, which is random 
enoug h on larger scales f or strong scattering to still be appli- 
cable dPogiel et all |1994|) . On the other hand, in some mod- 
els for ac celerating CRs beyond the knee (iLucek & Bell 120001: 
i Bell & L ucek. 2001), it has been argued that while downstream 
of a shock wave, strong turbulence is excited and isotropization 
of the CR distribution function is rapid enough to promote re- 
crossing, also no nlinear field amplification by CRs can occur up- 
stream (see also iMcKenzie & Volkl Il982) . Nonetheless, in our 
numerical calculations, we also have taken larger diffusion coef- 
ficients into account. As CR streaming generates waves for res- 
onant interaction, it has been suggested to assume that there is 
still a sufficient number of waves excited to guarantee pitch an- 
gle scattering, hence the validity of the first-order Fermi accel- 
eration process. As seen in the previous formula, the magnetic 
field strength B enters through the diffusion coefficient k(p), and 
according to V ■ B - 0, the adopted flux tube geometry (Eq. |7]l 
implies 



A{z)B{z) = AoBo, 



B{z) = Bo 



zo 



n-l 



(15) 



As the flow and the shocks in a flux tube are propagating along 
the magnetic field, we are dealing with parallel shocks here with 
no field compression. It has been suggested though that the max- 
imum energy gained by diffusive shock acceleration is higher for 
perpendicular shocks, when k is small, provided diffusion across 
the shock is fast en ough to allow for multiple shock crossings 
(e.g.'.TokiDii' Jl987h . 

It is straightforward to understand the acceleration charac- 
teristics, if we take a look at sufficiently large time scales f » f 
and therefore z s> zo- As seen in Fig. |6] the discontinuities are 
propagating at constant velocities, so that we can write, e.g., the 
distance of a shock from the galactic plane as z = u^t, where 
Ms denotes the shock velocity as given in Table [T] As a result, 
Eq. ( fT4l ) reduces via Eq. ( fT5l l for large z » Zo and p » 1 to 



K{p) 



1 mc 



3 ZeBozl 



u,t p . 



(16) 
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Fig. 7. The maximum momentum p^n^ in units of [mc] reached 
by a cosmic ray particle as a function of the shock distance Rs in 
[kpc] or time in [10*] years. Already within a few kpc from the 
galactic plane (corresponding to a flow time of about 10^ years) 
an almost constant particle momentum is achieved. The squares 
are calculated due to particle acceleration occurring only at the 
forward shock, and the diamonds also include the acceleration at 
the reverse shock. 

Approximating also Tacc - 4/c(/?)/Ms we end up with a simplified 
equation for the maximum momentum pm^x valid for the long 
term evolution 



dpn 



dt 



3ZeBozl 
Amc^ 



(17) 



This equation can be integrated easily, and specifying the initial 
conditions by p(fo) - Po at some time to, we get for t > to 



^^""^l^-7j' 



(18) 



which leads to an almost constant maximum momentum p^ax 
after the initial acceleration. 

This behavior of p^ax is clearly visible in Fig. Q where 
Eq. (fTSl l has been integrated in time for our velocity jumps. 
After the initial phase of a few 10'' years the momentum remains 
almost constant in accordance with the analytical estimate of 
Eq. (fTST l. As plotted in Figs. [3] and |5] the structure of the time- 
dependent flow exhibits two shock waves, and energetic particles 
can be (re)accelerated at both shock waves. The two curves de- 
pict the differences between the cases of one and two shocks. 
The full squares show the maximum momentum of particles if 
only the acceleration at the forward shock is considered; i.e., 
only the forward shock velocities ui and U2 are used in Eq. ( fT3l l. 
The upper curve reaching even a higher particle momentum is 
computed by using the velocity difference between the upstream 
region of the forward shock and the downstream region of the 
reverse shock. The acceleration at both shocks can therefore in- 
crease the maximum momentum by another factor of two com- 
pared to the acceleration only at the forward shock. From Fig.|2] 
it is evident that this reacceleration of energetic particles is very 
rapid and occurs early on within a few kpc above the galactic 
plane. 

Since particle acceleration here is essentially a first-order 
Fermi process, we expect the usual power-law spectrum to 
emerge. It is known from simple one-dimensional diffusive 



shock acceleration theory for steady plane shocks that for parti- 
cles injected at some momentum po, i.e. for a monoenergetic dis- 
tribution function, a power-law distribution for the downstream 
spectrum, fiip), will naturally arise, owing to the stochastic na- 
ture of the process, such that 



flip) 



No 3mi 

47r Ml — M2 



P \"l-"2 

po) 



(19) 



where A^o is the upstream particle number density with momen- 
tum p < Po- The upstream and downstream velocities mi_2 mea- 
sured in the shock frame are given by 



Ms 1 



Ml 



Ml 
U2 - : 

r 



(20) 
(21) 



with Mw the undisturbed wind speed, and r the compression ra- 
tio. The power-law index in Eq. ( fT9l l for a unmodified shock, 
i.e. in one where the back-reaction of the particles on the shock 
structure can be still neglected, is then 



? = 



3mi 



3(r, + 1)m2 



Ml - M2 2(M2 - 1) ' 



(22) 



where the upstream Mach number is given byM=^M-^j~ 
0.5—, for our standard case of zo - 15 kpc (see also Table [U; 
here Cj is the local speed of sound ahead of the shock front. 
The galactic wind is already highly supersonic near the disk, but 
the Mach number depends on the relative velocity between the 
shock and the wind, and only if the flow, which is catching up, 
is much faster than the quiescent wind, can we expect a strong 
gas shock and q x 4 during the first 20 million years. However, 
the compression ratio of the forward shock is r ~ 3 (s. Fig.|4]i, 
implying q - 4.5, which results in an energy spectral index of 
a ~ 2.9, including an energy dependent diffusion oc ^-O-^s Qf 
the particles propagating back to the disk. A mixture with CRs 
accelerated by the reverse shock with a lower compression ratio 
will tend to steepen the spectrum somewhat, which is already 
quite close to the observed value beyond the "knee". Since the 
injected particles are those coming from the disk, and already 
having a power-law spectrum, the resulting spectrum will also 
be a power-law. The spectral details will be the subject of fur- 
ther discussion in a forthcoming paper. 

5. Discussion and conclusions 

Soon after the diffusive shock acceleration theory was es- 
tablished as a majo r cand idate to expla in the origi n of cos- 
mic rays ('KrvmskV, 'TgTT'; ' Axford et all [19771 iBeUl ll978allH: 
Blandford & Ostriker, 1978), it became clear that supernova 
shock waves, which are the only non-exotic source of sufficient 
energy to provide the overall CR flux , woul d cut off in energy 
at about lO'"* eV (iLagage & Cesarskvl [T983I) . due to finite fife- 
time and curvature of the shock. A more recent analysis shows 
that energies up to Z ■ lO ' ^eV are possible, if diffusion is at 
the Bohm limit dBerezhkoi 11996 ). The observed spectrum be- 
tween the so-called knee and ankle thus have remained unex- 
plain ed for a long time. R ecent X-ray observations of Tycho's 
SNR dErikse n et al.Ll201 11) are interpreted such that particles are 
indeed accelerated up to the „knee" in regions of enhanced mag- 
netic turbulence corresponding to particles energies in the range 
of 10''* eV. However, the physical process of explaining energies 
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Table 1. Asymptotic wave velocities in [km/s] in the frame of the Galaxy and the maximal Mach numbers of the forward Mf r 
and reverse shock M,,,nax for time-dependent galactic winds, and the initial undisturbed wind velocity u„ at 1 Mpc. 



[kpc] 


[km/s] 


Vl 

[km/s] 


[km7s] 


Ik 
[km/s] 


11, 
[km/s] 


MfjTiax 
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10 


220 


710 


960 


880 


770 


80 


23 


15 


395 


780 


850 


900 


840 


47 


13 


20 


526 


840 


890 


980 


880 


46 
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of particles beyond the "knee" up to the "ankle" still remains 
unclear We have already mentioned the problems of detecting 
particles accelerated at the galactic wind termination shock. As 
discussed previously, several interesting suggestions have been 
put forward recently to explain the spectrum between the "an- 
kle" and the "knee . [Lu cek & Bell (2000) picked up on an ear- 
lier idea (e.g. JVolk & ]VIcKenziei.il98 1). arguing that the cosmic 
ray streaming itself would create sufficient turbulence upstream 
of the shock due to the rapid nonlinear growth of waves. If the 
SNR shock expanded i nto a stellar wind cav ity, a maximum en- 
ergy of ~ Z ■ 10"* eV (lBell&LucekLl200lh could be achieved, 
beyond which CRs are commonly believed to be of extragalactic 
origin. 



This hypothesis has received support in iDrurv et al.l (l2003b . 
They actually calculated the resulting spectra from a Sedov-type 
shock in the framework of a simple "box" model, in which the 
energetic particles are distributed uniformly within one diffusion 
length across the discontinuity, and are accelerated in momen- 
tum directly at the shock. They find a power-law spectrum with 
a small break at the knee all the way up to the ankle. The magni- 
tude of the break, however, depends on the diffusion length and 
on scales up- and downstream of the shock. It is therefore not 
well determined in such a simple model. The break occurs at 
the correct energy (where ordinary SNR shocks cease to accel- 
erate), but the power-law exponent in the momentum distribu- 
tion function can vary at the extreme between -3.5 (correspond- 
ing even to a hardening) and -5.75, while observations point to 
about -4.5. Therefore more detailed and presumably numerical 
calculations are needed in order to quantify the implications of 
this model. Re ferring to the self-consistent Galactic CR propaga- 
tion model bv iPtuskin etaP (ll997h [yolk & Zirakashvilil (|2004 
have questioned the suggestion by [Bell & Lucek (2001). They 
argue that the transition boundary between diffusion and advec- 
tion will have reached the termination shock at a distance of 
~ 300 kpc for particles with energies of ~ lO'^ZeV and that 
particles with higher rigidities diffuse freely upstream of the ter- 
mination shock and can only escape convectively downstream 
because of the small diffusion coefficient there, with no change 
in momentum. Thus the spectrum would be the same as the 
source spectrum in the disk, corresponding to a hardening of 
the ii'^-^-spectrum in the halo. In our model, acceleration oc- 
curs close to the disk, where the compression ratio is ~ 3, mak- 
ing the spectrum naturally steeper near the sources than in the 
disk. In addition, the particles have to diffuse back to the disk, 
leading to further steepening of the spectrum, as discussed in 
the previous section. A detailed discussion of this problem will 
be presented in a sequel paper Here we would just like to em- 
phasize that time-dependent changes in the boundary conditions 
of the galactic wind, naturally lead to the formation of strong 
shocks that will reacc elerate particles to en ergies well above the 
knee. We agree with IVolk & Zirakashvilil (l2004l) that the anal- 
ysis of the cosmic ray spectrum has to include the propagation 
effects of particles not only in the disk, but also in the Galactic 
halo. While we have analyzed these effects for starburst regions. 



it is clear that a similar scenario must hold for our Galaxy, al- 
though presumably restricted to smaller regions in the under- 
lying Galactic disk. Neverthele ss, 3D numerical simulations by 
lAvillez & Breitschw erdt (2004) have shown that outflow occurs 
even for thick Lockman and Reynolds type gas disks. The cru- 
cial point is that once overpressured hot regions have punched 
a hole in the lower halo, material can travel almost unimpeded 
to considerable heights, resulting in an outflow once the cosmic 
ray pressure gradient takes over These outflows will be modu- 
lated in a similar manner as the starburst types, due to changes 
in the star formation rate. The resulting shock waves can also be 
strong, because as we have emphasized, it is the relative velocity 
between fast and slow winds that matters, and steepening of the 
shock wave will occur as it runs down the galactic wind density 
gradient. 

In addition to the dynamical phenomena caused by SN- 
explosions, the magnetic field may also cause changes in 
the dynamics of the flow, and the idealized geometry will 
only approximate reality on average and on large scales. On 
smaller scales the field geometry can be quite complex, as 
shown by 3D numerical MHD simulations (Korpi et al., 199^ 
jAvillez & Breitschwerdtl l2005l e.g.). The tangling of field lines 
will inevita bly lead to a change in topology by magnetic re- 
connection (|Parkeri[l992h and to additional heating. Also MHD 
wave heating of the halo has been considered to be responsi- 
ble for highly ionized species (Hartquist, 1983). Field diffusivity 
and t urbulence, which are a nat ural consequence of a SN-driven 
ISM (lAvillez & BreitschwerdE 12005,) , will further promo te dy- 
namo action and B-field ampUfication dGressel et al.Ll2008l e.g.). 
Magnetized winds in other galaxies, e.g. dwarf galaxies, have 
also been modeled (Dubois & Teyssier, 2010) in order to study 
the enrichment of the intergalactic medium with magnetic fields 
within a cosmological context. 

Our model treats the base of the halo as an inner bound- 
ary condition, from which the flow emanates. In reality, there 
is a transition between disk and halo, and the outflow is a con- 
sequence of time-dependent SN activity in the disk. Again, a 
3D numerical simulation would be necessary to properly model 
fountain flows (Avillez & Breitschwerdt, 2004, e.g.). However, 
none of the 3D models put forward so far has considered the ef- 
fect of CRs and their nonlinear coupling to the flow, which is the 
purpose of this analysis. 

Although also neglecting effects of Galactic rotation and 
nonlinear wave damping, our model has the advantage that 
particle acceleration is calculated self-consistently along with 
the galactic wind flow. The hydrodynamic version of the trans- 
port equation that we use even allows us to enter the nonlinear 
regime, where the shocks can be modified by the increasing CR 
pressure, so we are not restricted to test particle arguments. Such 
a self-consistent picture allows us not only to deduce interesting 
quantities like mass-loss rate of a galaxy, but also to calculate 
e.g. X-ray spectra that can be compared to XMM-Newton 
and Chandra observations, in particular for delayed recombi- 
nation from nonequilibrium distributions of ions in the halo 
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(iBreitschwerdd. 120031 cf.)- In a next step we will include the CR 
momentum distribution function to numerically calculate the 
resulting cosmic ray spectra between the "knee" and the "ankle". 

To summarize, in this paper we have shown that 

1. time-dependent galactic winds have an asymptotic behav- 
ior that matches the st eady-state solutions obtained by 
iBreitschwerdt et al] (Il991h very well; 

2. these steady-state solutions are stable with respect to pertur- 
bations; 

3. changes in the inner boundary conditions are very likely to 
occur, since the flow time of the wind gas is on the order 
of 10^ years, so much longer than the switch-on time for 
superbubbles or starbursts; 

4. changes in the boundary conditions, in particular an increase 
in pressure due to increased star formation rate, result in 
time-dependent features like shocks (forward and reverse) 
and contact discontinuities; 

5. these shocks running down a density gradient can become 
very strong and propagate with roughly constant velocity; 

6. cosmic rays, propagating along with outflowing gas from 
the disk, can be accelerated very efficiently to energies up 
to 10'^ - lO'^ZeV by these shocks, thus suggesting a new 
mechanism for generating the Galactic cosmic rays between 
the "knee" and the "ankle"; 

7. particle acceleration is very fast and occurs close to the 
galactic disk, typically with a few kpc; 

8. the particle spectrum at the shock is most likely steeper than 
in the disk, because the Mach number of the forward shock 
depends on the relative velocity between the background 
wind and compressional waves, and is therefore smaller 
This leads to a spectral index < -4, i.e. a steepening at the 
shock. In particular, in the Milky Way, where bursts of star 
formation are less violent, shocks forming in the halo will 
have lower Mach numbers, thus leading to a steepening of a 
spectrum beyond the "knee"; 

9. in addition the spectrum must also steepen, because the parti- 
cles have to diff'use back to the disk; thus, the lower galactic 
halo seems to be a promising site for accelerating CRs be- 
yond the "knee". 
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Appendix A: Discretized equations 

The system of equations presented in Sect. 12.11 are rewritten to 
an integral form, which is necessary for applying a finite volume 
discretization. Following this procedure ensures the numerical 
conservation of mass, momentum, and energy. 

Introducing the temporal operator by & = jc - x"'"' between 
the two time levels t = f'"' + 6t, as well as the spatial difference 
Axi - xi+i - xi, we get a discrete version of the corresponding 
physical equations ([U - (|5]i. By defining the discrete volume bet- 
ween the distances zi and zi+\ through 



Az' 
AV, = Ao I Az, + -4 

JZfi 



we get, e.g.. 

dAV, ^ 
dzi+i 



i/+l 



An 1 + 



7^ ^ 



^0 ^ 



(A.1) 



(A.2) 



Since the adaptive grid changes the locations of grid points, we 
have to compute the relative velocity 

6t 



ui ■ 



(A.3) 



with respect to the grid motion and m'^''' appears in all advection 
terms. 

Written in a conservation form for a the flux-tube geometry, 
the discrete equation of continuity ([U reads as 

6 (piAVi) I5t + A [Aiufp) = , (A.4) 

where p; denotes the density upstream to the relative velocity 
M^'^'. Owing to the mathematical properties of hyperbolic con- 
servation laws, only those discretization schemes are stable that 
take th e direction of the flow into account(see e.g. LeVegue 
(1 19901) for all mathematical details). The equation of motion ^ 
is given by 

6 (piuiAVi) I St + A {Aiufpfii) + [P^j + P^j + Pw,/) AAi 

+ pi Ay;V<Dgai - MQ,;Ay, = , (A.5) 

where mq / stands for the artificial pressure term added to broaden 
shock fronts over a finite length scale, typically Iq - lO^^z. The 
gas energy density ^ reads in the volume-integrated discrete 
form 

6[E^jAv)l6t + A(A/Mf£g,,) + Pg,,A(A,M/) 

- eQ,,p,Ay, = , (A.6) 

which also includes the viscous energy dissipation term eq. The 
cosmic ray energy density (|4} is given through 

6{E,jAVi)l6t + A[AiufE,,) 

+ (rc-l)£c,/A(A,M,) = 0. (A.7) 



Finally, we have also to solve the equation of wave energy den- 
sity © 



5{E^jAVi)l6t + A[A,ufE^,^ 

+ (yw- 1 )£■«,; A(A/M/) = 0. 



(A.8) 



This set of physical equations is augmented by the grid equa- 
tion, which redistributes the grid points within the computational 
domain, and all details about this numerical technique can be 
found in lDorfi & Drurvl (11987 ). Therefore, in total we obtain 6A^ 
discrete nonlinear algebraic equations for the unknowns p/, m/, 
Eg,i, Ec,i, Ev,,i, and zi at every index I, where A^ = 300 is the 
number of spatial grid points. The resulting algebraic system is 
solved at every new time step by adopting a Newton-Raphson it- 
eration procedure where the corresponding Jacobian matrix has 
to be calculated from the discrete set of equations (IA.4t - (IA.8b . 

Due to the implicit nature of our numerical scheme we can 
compute time-dependent as well as stationary solutions, which 
can directly be compared to the stationary solutions obtained 
from solving the ODEs from the time-independent equations 
(Breitschwerdt et al., 1991). Furthermore, in the steady state 
case integrals of the equations can be derived and used for a 
numerical check of the code. Both methods show that the local, 
as well as global differences between the solution of the ODE 
system and the numerical scheme presented above are less than 
lO"*" everywhere. This small difference remains because of the 
staggered mesh representation of the physical quantities and the 
relative accuracy of the Newton-iteration of 10"^. 



Appendix B: Artificial viscosity 

Shock fronts are treated by appl ying an artificial viscosity in ten - 
sor form based on the work of iTscha rnuter & Winkler (1979!l. 
Since we are using an adaptive grid (jDorfi & Drurv. .1987 ) the 
length scale of the broadening /q can be specified according to 
our physical needs, i.e. /q <K Iqr = k/u^, which defines the length 
scale due to the particle acceleration. This way we can ensure 
that particles are accelerated by the 1. -order Fermi mechanism 
because the shock front is seen as a discontinuity by the ener- 
getic particles. 

The tensor formulation of the artificial viscosity requires 
some calculations for the adopted flux tube geometry since the 
pressure tensor Q, the viscous momentum transfer mq and the 
viscous energy d issipation sp have to be calcu lated. Without 
more details (see ITscharnuter & Winklen Il979b we state these 
quantities 



Q = ^Q 



(Vu)-iv-Ml 



MQ = -V-Q 

fiQ = --Q-(Vu), 
p 



(B.l) 
(B.2) 
(B.3) 



where (Vu) is the symmetrized velocity gradient and I the unity 
tensor 

Artificial viscosity coefficient fj.Q consists of a linear and a 
quadratic term weighted by qi and q2, i.e.. 



y^Q = -qilqCs + ?2^qmin(V-U'0) 



(B.4) 



The definition of /iq include a typical viscous length scale /q. 
When inspecting the last expression one sees that the linear term 
in Zq also scales with the sound velocity Cj and is always present 
when qi + 0. The second term is quadratic in /q, and it is evident 
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that compressive and nonhomologous motions; i.e., V-u < and 
fr(Q) ^ produce a viscous pressure. The amount of aitificial 
viscosity is then determined by these length scales q\l(^ and ^I'q- 
To obtain the corresponding expressions (IB. 2) and (IB. 3) for 
our flux tube geometry, we recall that the distance from the 
galactic plane is denoted by z and the area of the flux tube is 
defined by A(z). Neglecting the curvature within the flux tube 
area, we get the simple orthogonal metric tensor by 



'iik - 



( 1 











,2 





io 





A2 



(B.5) 



Following the rules derived bv JTschamuter & Wintderl (Il979l) . 
we can calculate the divergence of the velocity from the last ex- 
pression 

dAu 

as well as the symmetrized velocity gradient 
{ du 



1 dAu 
'^"A~d^ 



(B.6) 



^' = 



Tz ' 

u dA 
lAdl 

















u dA 

TA'dl 

The viscous pressure tensor is given through 
„ 1 dAu 



(B.7) 



A dz 
du 



X 

1 dAu 
3l"5F 






(B.8) 



u dA 
2A5F 



1 dAu 







u dA 
lAdl 



1 dAu 
3A~5F 



and has the desired property of being traceless; i.e., there is no ar- 
tificial viscosity in the case of homologous motions. According 
to the previous terms, the viscous pressure is given in the flux 
tube geometry through 



mq - 



1 d 

A3/2 5^ 



,1/2 



dAu I du 



^^'^\Tz-- 



1 dAu 
3~dr 



(B.9) 



The energy dissipation by the viscosity is determined from the 
contraction of the pressure tensor with the symmetrized velocity 
gradient (see Eq. IB.31 l. To ensure physical meaningful results it 
is necessary to find a discrete version, which always guarantees 
a positive viscous energy generation. This term has been imple- 
mented in our code through 



'^^-2- 



3 jUq dAu 
'dT 



du 

Tz 



1 dAu 
3A~5F 



(B.IO) 



guaranteeing that the viscous energy generation eg is also nu- 
merically positive. 

According to the discretization scheme, the artificial viscos- 
ity terms are written as 



P^QJ 



?i^qCs,;-^2/qmin — — ■ 



A(zfM/) 







(B.ll) 



including the linear part with Cj,/, the adiabatic sound veloc- 
ity at a grid point zi- The corresponding discrete and volume- 
integrated versions of mq and eg are then given as follows. The 



viscous pressure term in the discrete equation of motion 
reads as 



mq./AV/ 



(B.12) 






3/2 



m/P'- 



A(A;M;) 

AV, 



Am; 

Az, 



1 A(A/M/) 
3 



and the corresponding viscous energy dissipation 
cretized at a grid point zi according to 



AV, 

A.6t is dis- 



SQjpiAVi = -fiQjpiAHAiui)) 



Am/ 

Ai; 



1 AA;M/ 

3 AV, 



n2 



AVl. 



(B.13) 
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